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Abstract 
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1 Introduction 



There are a few processes measured at the LHC where the need for corrections 
beyond next-to-leading order (NLO) is free of doubt. One of them is top quark 
pair production, and the completion of a full NNLO prediction is well under 
way [T|2l3|^|5ll6ll7|8|9|T0|ll|12|[T3] . Soft gluon and Coulomb effects also have 
been taken into account beyond the next-to-leading logarithmic accuracy and 
have been combined with fixed order results to come up with predictions as 

Among the key ingredients of the full NNLO calculation are the master inte- 
grals entering the two-loop virtual corrections. While in [22|23f24] the latter 
enter in a semi-numerical form, a fully analytical representation for a subset of 
the needed master integrals has been presented in [HUE], where the fermionic 
and leading colour contributions to the qq channel and the leading colour con- 
tributions to the gg channel are calculated. An analytical representation for 
the non-planar seven-propagator integral occurring in the light fermionic cor- 
rection to the gg channel also has been achieved [7f8] . However, explicit results 
for some of the most complicated non-planar master integrals are still missing. 

Here we present numerical results for two non-planar seven-propagator topolo- 
gies, one entering the light fermionic correction to the gg channel, where ana- 
lytical results exist [7|8] , the other one entering the heavy fermionic correction 
to the gg channel, where no analytical result is available yet. Apart from scalar 
master integrals, we also give results for an irreducible tensor integral of rank 
two for the diagram entering the heavy fermionic correction to the gg channel. 
In order demonstrate the applicability of the tensor option to various types of 
integrals, we also calculate some two-loop two-point functions involving sev- 
eral different mass scales and show that the timings for the tensor integrals 
are not much larger than the ones for the scalar integrals. This means that 
a numerical approach in certain cases can help to alleviate or even avoid the 
procedure of amplitude reduction to master integrals. 

The results have been obtained with the program SecDec 2.1 [25|26f27j . 
based on sector decomposition to disentangle the singularity structure, fol- 
lowed by numerical contour integration. Compared to version 2.0 of SecDec , 
version 2.1 contains a number of new features, which are also presented in 
this article. Among them is the possibility to evaluate tensor integrals with 
(in principle) no limitation on the rank. Another new feature is the option 
to apply the sector decomposition algorithm and subsequent contour defor- 
mation on user-defined functions which do not necessarily have the form of 
standard loop integrals. In fact, to achieve a convenient representation for one 
of the seven-propagator integrals, some analytical manipulations have been 
done before starting the algorithm, this way reducing the number of produced 
subsectors, leading to improved numerical behaviour. 
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The structure of this article is as follows. In Section [2] we will derive the expres- 
sion serving as a starting point for the evaluation of the massive non-planar 
two-loop box diagram mentioned above, and describe novel types of transfor- 
mations which can be used to reduce the number of sector decompositions. 
The new features of the program SecDec 2.1 are discussed in Section [3) Nu- 
merical results are presented in Section HI where we also give results for some 
rank three tensor integrals for massive two-point functions. An appendix con- 
tains a manual-style description of the new features of the program. Detailed 
documentation of the program also comes with the code, which is available at 
http : / / secdec . hepf orge . org. 



2 Analytic preparation of the master diagrams 

The structure and usage of SecDec and the procedures it uses are described 
in detail in Refs. (25, 26 28 29]. The main purpose of this section is to explore 
the possibilities arising from a mixed approach, where simple analytic manip- 
ulations to the integral before feeding it to the sector decomposition algorithm 
can lead to a large gain in efficiency for the subsequent numerical evaluation. 

2. 1 Non-planar seven propagator integrals with two massive on-shell legs 



D B , 




(a) ggttl (b) ggtt2 

Fig. 1. Massive non-planar two-loop box diagrams entering the heavy (a) and light 
(b) fermionic correction to the gg channel; the thick lines denote massive particles. 

The diagrams shown in Fig. [1] are master topologies occurring in the two-loop 
corrections to ti production in the gg channel. While results for the diagram 
corresponding to massless fermionic corrections in a sub-loop - which we call 
ggtt2 - are available in terms of (9(800) generalized polylogarithms [7], ana- 
lytic results for the integral ggttl, containing a sub-diagram with a massive 
loop, are not available. Numerically however, the evaluation of ggttl is easier 
than the one of ggtt2, due to its less complicated infrared singularity structure. 
While the leading poles of ggtt2 are of order 1 / e 4 , and intermediate expressions 
during sector decomposition show (spurious) pole structures where the degree 
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of divergence is higher than logarithmic, the integral ggttl has only finite con- 
tributions. Therefore we evaluate ggttl with SecDec 2.1 fully automatically, 
while for ggtt2 it is advantageous to make some analytical manipulations be- 
forehand, reducing both, the number of Feynman parameters to be integrated 
over numerically and the degree of divergence. 

The expression for the scalar integral ggtt2 in momentum space is given by 

d D /ci d B k 2 



" ^i7T2^ J D Rl D R2 D Ri D Ri D Bl D B2 D B3 



where D = 4 — 2e. The Feynman propagators D Ri corresponding to the "rhom- 
bus" sub-loop in Fig. [I] are given by 



D Rl = (h - k 2 y + i5 , D R2 = (h -k 2 + pa) + iS 



D Rs = (k 2 + p 4 y + i5 , D R4 = (k 2 + pi + p 4 ) + iS , (2) 

where the p; are the external momenta with pg = p| = m? and p\ = p| = 0, 
and ki, k 2 are the loop momenta. Integrating out the loop momentum k 2 first, 
we are left with an expression containing k\ and external momenta only, to 
be combined with the propagators 

D Bl = (h - p 3 ) 2 + id , D B2 = (h + p 4 ) 2 + M , D B . A = k\ - m 2 + i8 . (3) 



The introduction of Feynman parameters for the one-loop subgraph X R con- 
taining only the loop momentum k 2 leads to 

= i / n n dD n n = ^ + c) / f[ ds,*(l - £ , 
17T i J D Rl JJ R2 L> R3 L> Ri J i=1 j=1 

(4) 

with 

—F{X, ki) = D Bl X x X 2 + [kx + Pi + PifxxXz + (fci + p 2 + P4) 2 ^2^4 + D B2 XzXi . 

We eliminate the 5-function in eq. (J4]) with the substitution 

Xt = t 2 (l - *a) , ar 2 = tit 3 , s 3 = (1 - ti)t 3 , (5) 

to achieve a factorisation of the parameter i 3 which then is integrated out 
analytically, leading to 

_ 2 r(2 + 6)r 2 (i-6) * fiikY^ 
Xr --1 — fxT^ — Jo dtl Jo dt ^^ 
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with 

-F(t, fa) = D Bl t x t 2 + (fa +pi+ p A ) 2 t x t 2 + (fa +P2+ P^) 2 tit2 + DB 2 iih , 

(6) 

and where we used the shorthand notation = 1 — tj. Now we combine 
the expression for the 1-loop rhombus X p with the remaining /^-dependent 
propagators, treating the expression of eq. flBJ) as a fourth propagator with 
power 2 + e, to obtain, after integrating out fa, 

2r(3 + 2e)r 2 (l -e) f\ f\ 
9 »'=1 r(l-fc) /o df '/o dfaX 

4 /•! 4 

-3-26 7/.._^U+3e 



4 ! 4 

]J / dzj ^4 +e 5(1 - J2 z i) Fmv(z, h, t 
i=i Jo ;=i 



2 )-^ Z<^(f) 1+ * , (7) 



where 



Utf-p(z) = ^^Zj and 
i=i 

Ttf V (z, ti) = -s u z 2 z 3 - - 5iZ 2 ^4 - S 2 z 3 z 4 + m 2 z 1 (z 1 + z A Q) , (8) 

with 

X = •§13^1^2 + S 23 tit 2 j Si = Si 2 tit 2 , S 2 = Si 2 tit 2 

Q = tj 2 + ti* 2 , % = (pi + pj) 2 . (9) 

Now we eliminate the 5-function by performing a primary sector decomposi- 
tion [28] in zi, . . . , Z4 to obtain 

2 r(3 + 26)r 2 (l-6) ri ri * , 
G NP =- ^^—^ J q dti j Q dt 2 £ G np , (10) 



=i 



with 



G l NP = [ 1 dz 2 dz 3 dz 4 z\ + ' (l + Z2 + z 3 + z 4 ) 1+3t F\Z,ti)- z -* 
Jo 

F l (z, ti) = —s\ 2 z 2 z 3 - Tz 4 - Siz 2 z A - S 2 z 3 z A + m 2 (l + z A Q) , 

G 2 NP = f dzidz 3 dz 4 zl +e (1 + zi + z 3 + z 4 ) 1+3e -F 2 (z,ti)~ 3 ~ 2e 
Jo 

JF 2 (z, ti) = -s l2 z 3 - Tz x z A - SiZ A - S 2 z 3 z 4 + m 2 zi(zi + z A Q) , (11) 

G 3 NP = [ d Zl dz 2 dz 4 zl +e (1 + Zl + z 2 + z 4 ) 1+3e T^U)- 3 - 26 
Jo 

J~ 3 (z, U) = -si 2 z 2 - Tz\z A - S\z 2 z A - 5*2-24 + m 2 zi(zi + z A Q) , 

G% P = [ 1 dz 1 dz 2 dz 3 (l + Zl + z 2 + z 3 ) 1+3e ^(z,U)- 3 - 2e 
Jo 

F 4 (z,ti) = -s 12 z 2 z 3 - Tzi - Siz 2 - S 2 z 3 + m 2 z x (z x + Q) . (12) 
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We observe that ^(zjU) is of the form m 2 + func(zj, £j), so does not need 
any further decomposition, and primary sector 3 can be remapped to primary 
sector 2 by exchanging z 2 -H- z 3 and Si <H> S 2 . Hence, we are left with the 
treatment of primary sectors 2 and 4 only. 

The integrals G^p can have singularities both at zero and one in t\ and t 2 . 
With the sector decomposition algorithm, only singularities at zero are fac- 
torized automatically. Consequently, we remap the singularities located at the 
upper integration limit to the origin of parameter space by splitting the inte- 
gration region at | and transforming the integration variables to remap the 
integration domain to the unit cube This procedure results in 12 integrals, 
some of which already being finite, such that no subsequent sector decompo- 
sition is required. Some of the integrals however lead to singularities of the 
type Jq 1 dx x~ 2 ~ € after sector decomposition, which we call linear singularities. 
These singularities are spurious and can be subtracted by expanding the Tay- 
lor series in the subtraction procedure up to the second term [2"gf3"U] . However, 
this procedure can lead to large cancellations between subtraction terms and 
therefore to numerical instabilities. Hence it is advisable to try avoiding this 
type of singularity from scratch. In the following section we describe a method 
which can help to do so. 

2.2 Removal of double linear divergences via backwards transformation 

The aim of the procedure described in this section is to achieve a transforma- 
tion of potential linear divergences into logarithmic divergences as far a pos- 
sible. A different procedure towards this goal, based on integration-by-parts 
identities, has been described in [25] . The latter method however can increase 
the number of functions to be integrated substantially, while the method de- 
scribed below in general reduces the number of further iterations and therefore 
the number of produced functions. Yet another method to reduce the number 
of functions produced during factorization has been suggested in Ref . [31] , but 
we do not use it here as it can introduce singularities at the upper integration 
limit which subsequently have to be remapped again. 

To explain the type of transformation advocated here, we use a function of 
the structure of T 2 in eq. ( ITT]) as an example, renaming t 2 —> z 2 ,ti — > z 5 . 
Concerning the Feynman parameters, we can identify the following structure 
in eq. ( ITT]) (in our concrete example N — 5, Zj — z±, z k — z\) 

I = IJ { f dzi) [zj (P(z jk ) + z k Q(z jk )) + R(z jk )}~ a , (13) 
i=i ^ Jo > 

where a is assumed to be positive, and where P, Q and R are polynomials of 
arbitrary degree of the Feynman parameters Zj k = (z\, z k , ... , z n) 
and kinematic invariants. The hat denotes Feynman parameters which do not 
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occur in the corresponding function. 

In eq. (fT3|) . all terms multiplied by the Feynman parameter z k are also multi- 
plied by the Feynman parameter Zj. Hence, the sector decomposition method 
can be applied "backwards". To explain this in more detail, consider the fol- 
lowing function: 

N f f 1 1 1 

J = II \J dz i \ — \ Z 3 P (Zjh) + z kQ{z jk ) + R{Zjk)] a [ ®(z k - Zj) + Q(Zj - Z k )} 

(14) 

Now we substitute Zj = z k tj in sector (1) and z k = Zjt k in sector (2), to 
obtain, after renaming again tj into Zi 

N r r 1 

J = II / dz i\- + ZkQ{z jk ) + R(z jk )]- a (15) 

t=i J Zj 

N ( r 1 11 

= II / d ^ - [zk{zjP{z jk ) + Q{z jk )) + i?(ij fc )]- Q (16) 

+ II { f [ z AP(z jk ) + z k Q(z jk )) + . (17) 

i=i uo J 

We observe that the term in eq. (fTTl) is the same as eq. (|T3l) . Therefore we can 
replace it by the expressions (|T5|) minus (|T6l) . The effect is twofold: The degree 
of the polynomial in ZjZ k is reduced in eq. fll5p . and in eq. f[T5]) the degree of 
divergence in Zj is reduced if a > 1. 

After all transformations of this type we arrive at a total of 15 functions partly 
needing an iterated sector decomposition. Together with the introduction of 
the new feature of user-defined functions in SecDec 2.1, which we will de- 
scribe in the following section, we are now able to compute the ggtt2 diagram 
in a reasonable amount of time. 



3 Structure and new features of SecDec version 2.1 



3. 1 Structure 



The workflow of the program is shown in Fig. [2J The directory structure of 
SecDec splits into two main branches: loop and general. In the loop part, 
multi-scale loop integrals can be evaluated without restricting the kinematics 
to the Euclidean region. Integrable singularities are dealt with by deforming 
the integration contour into the complex plane. In the general part, the 
poles of more general parametric functions can be factorized, and subsequently 
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loop directory 



Feynman loop integral 

(generated automatically) 



any integral matching 
loop integral structure 

(inserted by user) 



primary sector 
decomposition 



iterated sector 
decomposition 



multiscale? 



no 




iterated sector 
decomposition 




contour 
deformation 



multiscale? 



no 



subtraction of poles 



general directory 



more general 
parametric function 

(inserted by user) 



factorization 

optional: 

remapping of endpoint 
singularities 



I 



iterated sector 
decomposition 



expansion in £ 




numerical integration 




result: 






Laurent series in £ 



Fig. 2. Flowchart showing the main steps the program performs to produce the 
numerical result as a Laurent series in e. 

the functions can be evaluated numerically, provided they contain only end- 
point singularities. Contour deformation is not available in the subdirectory 
general. For more information on the structure of the various subdirectories 
we refer to Ref. 1261. 



3.2 New features 



The main new features of SecDec version 2.1 are the -u option in the loop 
directory to decompose "user denned" functions rather than standard multi- 
loop integrals, an extended tensor integral option, and improvements in the 
error treatment. 



Evaluation of user-defined functions with arbitrary kinematics 



If the user would like to calculate a "standard" loop integral, it is sufficient 
to specify the propagators, and the program will construct the integrand in 
terms of Feynman parameters automatically. However, in section [2] we have 
seen that an analytical step can be helpful when dealing with complicated inte- 
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grals. Integrating out one Feynman parameter analytically reduces the number 
of integration variables for the subsequent Monte Carlo integration and there- 
fore in general improves numerical efficiency This implies that the constraint 
5(1 — J2i x i) has been used already to achieve a convenient parametrisation, 
and therefore no primary sector decomposition to eliminate the 5-constraint 
is needed anymore. In such a case, the user can skip the primary sector de- 
composition step and insert the functions to be factorized directly into the 
Mathematica input file, using his favourite parametrisation. More generally 
speaking, the purpose of this option is to be more flexible with regards to 
the functions to be integrated, such that expressions for loop integrals which 
are not in the "standard form" - for example due to analytic manipulations 
which have been performed already on the integral - can be dealt with as 
well. This includes the possibility to perform a deformation of the integration 
contour into the complex plane, taking the user-defined functions as a starting 
point. Oriented at the functions T and U for the "standard" loop case (see e.g. 
eq. ©), the user-defined functions can encompass the product of two arbitrary 
polynomial functions with different exponents and an additional numerator. 
Details about the usage of this option are given in appendix IA.2I 

The tensor integral option and the syntax for the definition of the numerator 
in terms of contracted loop momenta are described in detail in appendix IA.3I 



Error estimates 

When dealing with complicated integrands it can happen that the error given 
by the numerical integration program - which is based on the number of sam- 
pling points only - underestimates the true error. The numerical integrators 
contained in the Cuba library [32|33] give an estimate of the correctness of 
the stated error to a given result. The new SecDec version 2.1 collects the 
maximal error probability for each computed order in the dimensional regula- 
tor e and writes it to the result files * . res. In the generic case, the information 
on the reliability of the stated error is given as a probability with values be- 
tween and 1. If the integrator returns a value larger than one, the integration 
has not come to successful completion, and a warning is written to the result 
files. This feature helps the user to judge how reliable the numerical results, 
respectively the given errors, are. 



3.3 Installation and usage 

The program can be downloaded from http : / / secdec . hepf orge . org. 
Unpacking the tar archive via tar xzvf SecDec.tar.gz will create a directory 
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called SecDec with the subdirectories as described above. Installation is done 
by changing to the SecDec directory and running ./install. 

Prerequisites are Mathematica, version 6 or above, Perl (installed by default 
on most Unix/Linux systems), a C++ compiler, and a Fortran compiler if the 
Fortran option is used. 

Details about the usage can be found in the appendix, where manual-type 
descriptions of the new features are given, and also in Ref. [26] and the docu- 
mentation coming with the program. 



4 Results 



Now we present numerical results for two non-planar seven propagator master 
topologies occurring in the two-loop corrections to the process gg — > ti, shown 
in Fig. [TJ In addition to the scalar integrals, we also compute an irreducible 
tensor integral. In order to demonstrate the applicability of the tensor option 
to various contexts, we also give results for some rank three tensor two-point 
functions involving several mass scales. 



4-1 The ggttl diagram 



Numerical results for the diagram ggttl (see Fig. Ula)) are shown in Fig. |3]for 
both the scalar integral and an irreducible rank two tensor integral. 



--++++++ 



Non-planar scalar massive 2L box diagram ggttl 

* s» — ** 



Real (SecDec] 
Imag (SecDec] 



(a) scalar integral 



Non-planar rank 2 massive 2L box diagram ggttl 



Real (SecDec! 
Imag (SecDec) 



(b) rank 2 tensor integral 



Fig. 3. Results for the scalar integral ggttl shown in Fig.[T|a), and the corresponding 
rank two tensor integral ggttl with k\ ■ hi in the numerator. We vary s\2 and fix 
•S23 = -1-25, mi = m\,pl = p\ = m\ = 1. 
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The integral representation of the diagram ggttl is given by 

Gggtn = (t4) 2 / d n hdD n ^ D ^ = k l~ < D 2 = (h + Pl ) 2 - m\, (18) 

B z = k\- m|, D 4 = (k 2 + p 2 ) 2 - ml, D 5 = (ki - k 2 + pi) 2 , 
D 6 = (h -k 2 - p 2 ) 2 , D 7 = (fcx - k 2 + pi + p 3 ) 2 - ml , 

where we omitted the infinitesimal i5 in the propagators, and we use the 
convention that all external momenta are ingoing. The momenta ^3 and p± 
are massive on-shell: p\= p\ = m\. The results shown Fig. [3]^b) correspond to 
the rank two tensor integral with the same propagators and a factor of k\ ■ k 2 
in the numerator. The numerical integration errors are shown as horizontal 
markers on the vertical lines. The absence of such markers means that the 
numerical errors are smaller than visible in the plot. 

The timings for one kinematic point for the scalar integral in Fig.[3^a) range 
from 11-60 sees for points far from threshold to 1.6 x 10 3 seconds for a point 
very close to threshold, with an average of about 500 sees for points in the 
vicinity of the threshold. A relative accuracy of 10 -3 has been required for the 
numerical integration, while the absolute accuracy has been set to 10 -5 . For 
the tensor integral, the timings are better than in the scalar the nu- 

merator function present in this case smoothes out the singularity structure. 
A phase space point far from threshold takes about 5-10 sees, while points 
very close to threshold do not exceed 3600 sees for the rank 2 tensor inte- 
gral, Fig.^b). The timings were obtained on a single machine using Intel i7 
processors and 8 cores. 

For the results shown in Fig. Owe used the numerical values m\ = m\ = m 2 = 
1; S23 = —1.25, S13 = 2 m 2 — S12 — S23. We set m\ = m 2 for the results shown 
in Fig. [3] because this is the only case occurring in the process gg — > ti at two 
loops if the 6-quarks are assumed to be massless. Some numerical results for 
kinematic points with mi 7^ m 2 are shown in Fig. H] in order to demonstrate 
that it is easy to add another mass scale in our approach, while this would 
complicate analytical calculations enormously. 

4-2 The ggtt2 diagram 

Analytic results for the pole coefficients of the 1/e 4 and 1 /e 3 part of the dia- 
gram ggtt2 shown in Fig. [TJb) have been given in [7j. We confirm these results 
and show a comparison between analytic and numerical results in Fig. [5J 
The integral representation of the diagram ggtt2 is given by eq. ([TBI with 
m 2 = 0. For the results shown in Figs. [5] to [7| we used the numerical values 
p\ = p\ = m 2 = 1,S23 = —1-25, S13 = 2m 2 — s\ 2 — S23, and we extract an 
overall factor of 16 T(l + e) 2 . 
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Non-planar massive 2L box diagram ggttl , m 1 =1 



Real (SecDec) i 
Imag (SecDec) i 



Fig. 4. Results for the scalar integral ggttl shown in Fig. []Ja), with two different 
masses. We vary rri2 and fix s\2 = 5, S23 = — 1.25,p| = p\ = m\ = 1. 



Non-planar scalar massive 2L box 



Real (Analytic] 
Imag (Analytic) 
Real (SecDec) ^ 
Imag (SecDec) t 



4*\ 



Non-planar scalar massive 2L box diagram ggtl2 



.-+-+-++-+-+- 



Real (Analytic) ■ 

Imag (Analytic) ■ 

Real (SecDec) e 

Imag (SecDec) h- 



(a) leading pole 



(b) subleading pole 



Fig. 5. Comparison of (a) the leading and (b) the subleading pole coefficients be- 
tween the analytic result from [7] and the SecDec result, for both real part (blue) 
and imaginary part (red). We vary s\2 and fix S23 = — 1.25,p| = p\ = m 2 = 1- 

For two selected phase space points, we also compared the full result for all 
Laurent coefficients including the finite part to the result obtained analyt- 
ically [31] and find agreement within the numerical integration errors. Our 
numerical results for the remaining pole coefficients and the finite part of the 
diagram ggtt2 are shown in Figs. [6] and [71 

The timings for the leading and subleading pole coefficients of the diagram 
ggtt2 are ranging between fractions of a second and about 20 seconds. The 
coefficients of the 1/e 2 pole take between 13 and 300 seconds, while the coef- 
ficients of the 1/e pole take between 75 and 3000 seconds, depending on their 
distance to thresholds. For the finite part, the integration times range from 
250 seconds to 4000 seconds. For all Laurent coefficients, a relative accuracy 
of 5 x 10~ 3 has been required, which has not always been reached for the 1/e 
and e° coefficients. It also should be noted that the timings for points close 
to threshold depend rather sensitively on the settings for the Monte Carlo 
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integration parameters. 



f + 


Real (SecDec) I [— I 

Imag (SecDec) i 1 1 


- V + 
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^ ++¥ + + + + + + - 


- + + + 




- 1 + 
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+ T + 

R 



Real (SecDec) t 
Imag (SecDec) i- 



#N* +4+ + + + + + + 



(a) 1/e 2 pole 



(b) 1/e pole 



Fig. 6. Results for (a) the 1/e 2 and (b) the 1/e coefficients of the integral ggtt2. The 
kinematics are the same as in Fig. [5j 
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Real (SecDec) I 1 
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Imag (SecDec) i 1 1 
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Real (SecDec) I 
Imag (SecDec) h 



(a) finite part, including threshold 



(b) finite part, beyond threshold 



Fig. 7. Results for the finite part of the scalar integral ggtt2, (a) for a larger kinematic 
range, (b) zoom into a region further away from threshold. The kinematics are the 
same as in Fig. [5j The vertical bars denote the numerical integration errors. 



4-3 Massive tensor two-loop two-point functions 



Here we show that the option to evaluate integrals with a non-trivial numer- 
ator can also be applied to calculate two-loop two-point functions involving 
different mass scales, without the need for a reduction to master integrals. This 
fact can be used for instance to calculate two-loop corrections to mass param- 
eters in a straightforward way. As an example we pick the diagram shown in 
Fig. El 

We calculate scalar integrals and rank three tensor integrals, for two cases, 
one where m 3 is zero, and one where m 3 is nonzero. The tensor integral is 
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Fig. 8. Two-loop bubble diagram with different masses 



given by 



= {^) J D^D 5 5 (19) 

Di — k\ — m\, D 2 = (h + Pl ) 2 - m\, D 3 = (k, - k 2 ) 2 - m\, 

Da = (k 2 + pi) 2 - ml, D^ — k\ — m\. 



The fact that this tensor integral is reducible does not play a role here, because 
our purpose is to demonstrate that reduction may become obsolete considering 
the short integration times for the tensors. Results for the scalar and tensor 
integrals with m 3 = are shown in Fig. [HI while results for m\ = 3 are shown 
in Fig. dOj 



Scalar 2L bubble with m^=2.0, rri2=4.0 Rank 3 2L bubble with m^=2.0, rri2=4.0 




(a) scalar integral (b) rank 3 tensor integral 

Fig. 9. Results for real (blue) and imaginary (red) parts of the rank 3 two-loop 
bubble diagram shown in Fig. [HI (a) scalar case, (b) with numerator (k\ -k 2 ) {k\ -pi). 
The masses are m 2 = 2, m 2 = 4, m-s = 0. 

The timings for the massive two-point integrals are shown in Figs. [H](a) and 
(b), for m\ = and m 2 = 3, respectively. The timings were obtained on com- 
puters with Intel i7 processors and 8 cores. In both cases, a relative accuracy of 
0.1% was required for the Monte Carlo integration. In Fig.[TT](a), an absolute 
desired accuracy of 10~ 8 was set, in Fig.[TT](b) an absolute accuracy of 10 -6 . 
Fig. rm shows that for values of p 2 above the mass threshold at p 2 = ^.m\, the 
timings for the contracted rank three tensor integrals do not differ much from 
the ones for the scalar integrals. In the region below threshold, the timings are 
higher because the imaginary part is zero, and a vanishing function is difficult 
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to integrate numerically. As the relative error to a zero value is always infinite, 
the numerical integrator in this case tries to reach the absolute accuracy goal. 
If in addition to the vanishing imaginary part the real part is also close to zero, 
the integration times are highest, as can be seen from Figs. [TO] and [TT](b). To 
avoid artificially large timings in kinematic regions where the imaginary part 
is known to be zero, one could set contourdef =f alse for the kinematic points 
below threshold. This would have the effect that the imaginary part is not cal- 
culated at all, but set to zero from scratch, thereby reducing the integration 
time considerably. The minimal numerical integration time for a kinematic 
point above threshold in the case of the scalar two-loop bubble integrals is 
0.03 sees for m§ 7^ and 0.02 sees for m| = 0. 



Scalar 2L bubble witb m 2 =2.0. rri2=4.0, rri3=3.0 Rajik 3 2L bubble with m 2 =2.0, rri2=4.0. 013=3.0 




5 10 15 20 5 10 15 20 

P ! P 2 

(a) scalar integral (b) rank 3 tensor integral 



Fig. 10. Results for the rank 3 two-loop bubble diagram with three non-van- 
ishing masses, shown in Fig. El (a) scalar case, (b) tensor case with numerator 
(k\ ■ fcj) {k\ ■ pi). The masses are m\ = 2, m| = 4, m| = 3. 



Timings tensor/scalar 2L bubble with m,=2.0 r rri2=4.0 



(a) m 3 = 



Timings tensor/scalar 2L bubble with m 1 =2.0, m 2 =4.0, m 3 =3.0 



++ 

+++ 



(b) m 



Fig. 11. Comparison of evaluation times between the scalar and the rank three 
tensor integrals corresponding to a two-loop two-point function with 3 different 
masses. The red points are the evaluation times in seconds for the scalar integral 
at a given kinematic point, the timings for the rank 3 tensor integral are marked in 
blue. 
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5 Conclusions 

We have presented numerical results for two massive non-planar seven prop- 
agator topologies, one entering the light fermionic two-loop correction to the 
gg — > ti channel, the other one entering the heavy fermionic correction to this 
channel. For the latter, no analytical result is available yet. Apart from the 
scalar master integral, we also give results for an irreducible tensor numerator 
of rank two for this diagram. Our numerical results have been obtained with 
the program SecDec2.1, which is publicly available at 

http : / / secdec . hepf orge . org| Compared to version 2.0 of SecDec, version 
2.1 contains a number of new features, among them the possibility to evaluate 
tensor integrals with no principle limitation on the rank. The applicability 
of the tensor option to various types of integrals is further demonstrated by 
a number of results for two-loop two-point functions involving several differ- 
ent mass scales, where no analytical results exist. Another new feature is the 
option to apply the sector decomposition algorithm and subsequent contour 
deformation on user-defined functions which do not necessarily have the form 
of standard loop integrals. This new feature is used in combination with a 
novel type of analytic transformations, which can serve to reduce the number 
of functions to be integrated numerically. We believe that SecDec version 
2.1 brings us a major step forward in moving from the calculation of mas- 
ter integrals to the calculation of two-loop corrections for phenomenological 
applications. 
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A Appendix 

A.l Usage of the program 

(1) Change to the subdirectory loop or general. The setup in the loop 
directory should be used to calculate loop integrals and integrals with 
a similar structure. The option to use contour deformation is available 
for all functions processed within the loop directory. The setup in the 
general directory allows to evaluate more general parameter integrals, 
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which can have endpoint singularities at zero and one, but does not offer 
contour deformation. 

(2) In the general directory, copy the files param. input and 
template. m to create your own parameter and template files 
myparamf ile . input, mytemplatef ile .m and edit them according to the 
function to be calculated. In the loop directory, edit the files 
paramloop . input and templateloop .m if you want to compute a Feyn- 
man loop integral in a fully automated way. If you would like to de- 
fine a set of own functions rather than a standard loop integral, use the 
files paramuserdef ined. input and templateuserdef ined.m as a start- 
ing point. 

(3) Set the desired parameters in myparamf ile . input and specify the inte- 
grand in mytemplatef ile .m. 

(4) Execute the command ./launch -p myparamfile. input -t mytemplatefile.m 
in the shell. If you add the option -u in the loop directory, user defined 
functions are computed. 

If you omit the option -p myparamfile. input, the file param. input will be 
taken as default. Likewise, if you omit the option -t mytemplatefile.m, the 
file template. m will be taken as default. If your files myparamfile. input, 
mytemplatefile.m are in a different directory, say, myworkingdir, use the 
option -d myworkingdir, i.e. the full command then looks like ./launch 
-d myworkingdir -p myparamfile. input -t mytemplatefile.m, executed from 
the directory SecDec/loop or SecDec/general. 

(5) Collect the results. Depending on whether you have used a single machine 
or submitted the jobs to a cluster, the following actions will be performed: 

• If the calculations are done sequentially on a single machine, the results 
will be collected automatically (via the corresponding results*. pi 
called by launch). The output file will be displayed with your spec- 
ified text editor. 

• If the jobs have been submitted to a cluster, when all jobs have finished, 
execute the command ./results. pi [-d myworkingdir -p myparamfile] in 
the general, and ./resultsloop.pl [-d myworkingdir -p myparamfile] or 
./ resultsuserdefined.pl [-d myworkingdir -p myparamfile] in the loop di- 
rectory, respectively. This will create the files containing the final results 
in the graph subdirectory specified in the input file. 

(6) After the calculation and the collection of the results is completed, you 
can use the shell command ./launchclean [graph] to remove obsolete files. 

A. 2 Evaluation of user-defined functions in the loop directory 

In the following, we will describe the input and syntax needed for the files 
mytemplatefile.m and myparamf ile . input when treating functions which 
are different from the standard Feynman parameter representation which - in 
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the default setup - is derived automatically from the propagators or vertices 
of a multi-loop integral. 

The file mytemplatef ile .m should contain the following information: 

• List of user-defined functions: 

The user-defined functions should be polynomial in the Feynman parame- 
ters. They can contain monomial factors of Feynman parameters with arbi- 
trary exponents, functions of type U and T (i.e. polynomials in the Feynman 
parameters involving also kinematic invariants) with arbitrary exponents 
and a "numerator" with positive exponents only. An iterated sector decom- 
position is applied to U and T if the decomposition flag (see below) is set to 
B. In case the functions U and J 7 contain thresholds, a deformation of the 
integration contour into the complex plane becomes necessary. The setup is 
such that the integration contour will be formed based on the function J 7 . 
The list of user-defined functions must be inserted into mytemplatef ile .m 
using the following syntax 

functionlist= {function_ 1 Junctiori-2, . . . Junction J, ...}; 
with 

functionJ={# of functionalist of exponents} , {{function U , exponent ofU, 
decomposition flag} , {function J 7 , exponent of T ', decomposition flag}}, nu- 
merator}. 

The # of function of each user-defined function is a label of the function by 
an integer, where the default is just sequential numbering. However, there 
is also the option to label a set of functions with the same integer. In this 
case the functions are decomposed individually, but will be combined after 
the decomposition. This leads to fewer or simpler integrand functions when 
symmetries are found within a sector. It should be noted however that only 
functions which share the same exponent for all functions of type 1A re- 
spectively J 7 can be grouped together and therefore have the same function 
label. 

Each entry in the comma separated list of exponents corresponds to an 
exponent of a Feynman parameter occurring as a monomial in the Feynman 
integral. The decomposition flag should be A if no iterated sector decompo- 
sition is desired, and B if the function needs further decomposition. The 
numerator may contain several functions, with different non-negative expo- 
nents. 

• Dimension, kinematic conditions: 

The space-time dimension should be specified in mytemplatef ile .m via 

Dim= dimension. The default is Dim=4 — 2e. 

On-shell conditions can be specified in the list onshell={}. 

• Computation of the exponents of T and U (optional): 

If the exponents of the functions J 7 and U should be computed by the 
program according to the rules for standard multi-loop integrals, the user 
needs to specify the number of propagators and their powers in a list 
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powerlist=Table [power, {i , #propagators}] ; and the tensor rank of the 
diagram via ra.nk=rank; . 

Here, ^propagators should correspond to the number of propagators of the 
original Feynman diagram, or, more general, to the number of integration 
variables plus one. 

In addition to the changes made to mytemplatef ile .m, it is possible to set a 
maximal number of Feynman parameters occurring in the user defined func- 
tions by initalizing feynpars=. . . in myparamf ile. input. The default is the 
number of propagators N subtracted by one. 

A full example including detailed comments comes with the download of the 
program and can be found in /loop/templateuserdef ined .m and 
/loop/paramuserdef ined. input. 

A. 3 Evaluation of tensor integrals 

For the computation of tensor integrals, where the tensor is contracted with 
external momenta and/or loop momenta, the construction of the Feynman 
integral via topological cuts needs to be switched off, which corresponds to 
cutconstruct=0. 

Regarding the user input, the only additional information needed is the numer- 
ator in mytemplatef ile .m. Each scalar product of loop momenta contracted 
with either external momenta or other loop momenta should be given as an 
entry of a list: 

numerator={]?re/ 'actor , comma separated list of scalar products] . For example, 
a numerator of the form -2 k\-k 2 , where k\ and k 2 are loop momenta, should be 
given as numerator={— 2, kl * k2 }. A numerator of the form 2 (&i -pi){ki ■ k 2 ), 
where p\ is an external momentum, should be given as 
numerator={2, kl * pi, kl * k2 }. 
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